Genetic embedded matching approach to ground states in continuous-spin systems 
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Due to an extremely rugged structure of the free energy landscape, the determination of spin- 
glass ground states is among the hardest known optimization problems, found to be ./VP-hard in the 
most general case. Owing to the specific structure of local (free) energy minima, general-purpose 
optimization strategies perform relatively poorly on these problems, and a number of specially 
tailored optimization techniques have been developed in particular for the Ising spin glass and 
similar discrete systems. Here, an efficient optimization heuristic for the much less discussed case 
of continuous spins is introduced, based on the combination of an embedding of Ising spins into the 
continuous rotators and an appropriate variant of a genetic algorithm. Statistical techniques for 
insuring high reliability in finding (numerically) exact ground states are discussed, and the method 
is benchmarked against the simulated annealing approach. 

PACS numbers: 75.50.Lk, 02.60.Pn, 75.10.Hk 
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I. INTRODUCTION 



Complex (free) energy landscapes featuring a multi- 
tude of local minima separated by energy barriers are 
common in problems of statistical mechanics, chemical 
and biophysics consequently often subsumed under the 
label of "complex systems", be it biopolymers, struc- 
tural or spin glasses [1]. The consequences of this de- 
viation from the classical textbook situation of a poten- 
tial energy with at most a handful of metastable states 
are dramatic for the static and dynamic behavior of the 
affected systems in nature, including for instance "mem- 
ory" and "rejuvenation" effects in spin glasses Q, but no 
less pronounced for the theoretical investigation of mod- 
els of such situations with computational simulation or 
optimization methods: here, model systems get trapped 
in local minima for exponentially long times, prevent- 
ing an equilibration in finite-temperature simulations Q 
or lead to a vastly increased effort needed for an op- 
timization procedure to yield ground states with finite 
probability [H . Clearly, the presence of many local min- 
ima alone is not sufficient to pose serious problems for 
any optimization method more elaborate than a strictly 
downhill, local search. Rather, it is the organization of 
minima and interjacent barriers that is the cause for the 
trapping phenomena, and distinguishes the milder from 
the more severe cases While, for instance, many typ- 
ical biopolymers exhibit landscapes with moderate barri- 
ers separating minima of substantially different energies, 
with a "funneling" towards a unique global minimum Q , 
disordered and frustrated magnetic systems are rather 
characterized by many (quasi-) degenerate minima close 
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to the ground state(s) separated by large barriers, lead- 
ing to much more severe effects of metastability and slow 
relaxation @,@|. 

Independent of this connection between the structure 
of the energy landscape and the real or artificial (Monte 
Carlo or optimization) dynamics of a frustrated system, 
the problem of finding ground states or, alternatively, 
partition functions, of the corresponding models has been 
considered as a problem in the field of computational 
complexity Q. In computer science, traditionally mostly 
the "worst case" complexities of algorithms have been 
considered, i.e., the asymptotic scaling of the run time 
T(N) with the problem size N for the "hardest" set of in- 
put data within the class of allowed inputs, T max (iV) [f|. 
Quite generally problems with an asymptotically polyno- 
mial form of T' max (iV) are considered tractable, whereas 
an exponential divergence for the best known algorithm is 
associated with intractability. Paradigmatic results have 
been found for decision problems with "yes" or "no" an- 
swers, for which a powerful classification scheme has been 
established: problems with a known polynomial algo- 
rithm are grouped in V, whereas a more general class 
of problems for which the correctness of a solution can 
be checked in polynomial time is denoted MV. The po- 
tential hardness of the latter must then exclusively re- 
sult from the exponential growth of the search space, 
such that a theoretical computer capable of infinite par- 
allelism can solve such problems in polynomial time 0]. 
The hardest MV problems, namely those whose poly- 
nomial solution would imply polynomial complexity for 
all other MV problems, are termed MV complete, which 
includes most of the well-known hard problems such as 
the traveling salesman problem or the satisfiability prob- 
lem. While it is possible that all such problems might 
have polynomial-time solutions, i.e., V = MV, this is 
now considered to be extremely unlikely, and MV prob- 
lems almost certainly require an exponential computa- 
tional effort [8j. This classification extends to optimiza- 
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tion (instead of decision) problems where, specifically, 
those where the problem of deciding about the existence 
of a solution better than a given bound in the cost func- 
tion is MV complete, are termed MV hard. 
For the Ising spin glass with Hamiltonian Q 



H 



Si • Sj, 



(i) 



the problems of computing ground states or the partition 
function are known to be MV hard in space dimensions 
d > 3 or for two-dimensional (2D) systems in a magnetic 
field [h], [H| • The zero-field 2D problem, on the other 
hand, is tractable in polynomial time [H], [IH, [HI- In 
particular, ground states on planar graphs can be found 
by means of the mapping to a minimum-weight perfect 
matching problem, as discussed below in Sec. HU While 
thus, generically, spin-glass ground-state problems are 
hard, one has to keep in mind that this classification con- 
cerns the worst-case behavior among all possible realiza- 
tions of couplings Jij of the chosen distribution (e.g., bi- 
modal or Gaussian), whereas it is, for instance, simple to 
specify the ground state of the purely ferromagnetic sys- 
tem with Jij = Jo > 0, which also belongs to the allowed 
Jij realizations. Hence, relevant for actual computations 
is also the average complexity, depending on a chosen 
probability distribution P({Jij}) of couplings. Within 
the spin-glass phase, however, also this mean complexity 
is exponential for known exact approaches to the prob- 
lem in d > 2 [la ]. Correspondingly, heuristic optimiza- 
tion techniques for finding low-lying or ground states are 
called for. These might include generic approaches, such 
as simulated annealing [j^ . multicanonical [13] or par- 
allel tempering jl8} Monte Carlo simulations, but also a 
number of methods specifically tailored to the problem 
0, US HH, l2 3 . [23l. the latter generally showing the best 
performance [3|, [24| . Insofar as these methods make use 
of some type of relaxational (quasi-) dynamics, they to 
some extent also suffer from the slowness of relaxation 
entailed by the structure of frcc-cncrgy minima and sep- 
arating barriers. It should be pointed out, however, that 
such slow dynamics is not equivalent to hardness in the 
classifications of computational complexity [25| . Instead, 
slower than power-law relaxation of local dynamics also 
occurs in computationally polynomial systems J5| , such 
as the Ising spin-glass model in two dimensions [121 ] . The 
stochastic nature of most of these approaches requires 
a different description of their time complexity or effi- 
ciency: since such methods do not guarantee to yield 
ground states, one should now rather ask for the worst 
case or mean computational effort to end up in a ground 
state with an a priori prescribed success probability p s 
(for p s = 0.95, say), or for the distribution (over disorder 
realizations) of such minimal running times at fixed p s . A 
framework for such considerations will be developed be- 
low in Sec. IIVI Ground-state searches for spin-glass sys- 
tems are additionally complicated by an extraordinarily 
broad distribution of "hardness" over disorder samples, 
which draws into question the treatment of all samples 



with constant computational effort. In this context, it is 
discussed below in how far properties of individual disor- 
der samples can serve as hardness indicators and hence 
an automatic effort adaptation can be achieved. 

Ising spin-glass ground states have been considered 
with the aim to understand the nature of the low- 
temperature phase while avoiding the equilibration prob- 
lems of finite-temperature simulations. Ground-state 
computations for systems with different boundary condi- 
tions or with some fixed spins allow for the direct inves- 
tigation of domain- wall and droplet defects, whose prop- 
erties should reveal in how far finite-dimensional spin 
glasses are correctly described by mean-field theory (see 
Ref. [26| for a review of recent developments) . The poly- 
nomially tractable 2D case, in particular, has provided 
a fruitful playground for testing theoretical pictures of 
the spin-glass phase, and remains a topic of active re- 
search to date [27], HI, H^, [3(J. In terms of spin-glass 
phases realized experimentally, in particular in the mul- 
titude of systems with frustrating lattice structures that 
have come into focus more recently 3l|, systems with 
continuous spins are probably more common than the 
extremely anisotropic Ising case. In computing ground 
states for such systems, modeled, say, by the Edwards- 
Anderson Hamiltonian ([T]) with continuous O(n) spins 
Si, one leaves the relatively well- understood field of com- 
binatorial optimization. To my understanding, nothing 
is known about the (suitably generalized) computational 
complexity of this problem. It is easily seen [3, 8], how- 
ever, that already the g-state Potts spin glass corresponds 
to a multi-terminal flow problem known to be MV hard 
even in two dimensions for q = 3 [32j]. Nothing would 
seem to indicate that the XY case of continuous planar 
spins, or the Heisenberg model of 0(3) rotators could be 
easier computationally than the discrete Potts approxi- 
mation. With the exception of a study of the XY spin 
glass in the Coulomb gas representation [3^], all studies 
of low- lying metastable states in O(n) spin glasses (with 
n > 1) have relied on variants of a simple spin-quench 
technique corresponding to a T = Monte Carlo simula- 
tion with local updates [33,[35|,[36|,[37j (apart from studies 
of the computationally simpler case of the n — > 00 spher- 
ical spin glass [38| ) . This spin quench follows from noting 
that a necessary condition for metastability is that each 
spin be parallel to its local molecular field, 



hj = j, 



1 j Sj , 



(2) 



leading to the prescription of an iterative alignment of 
single spins Si parallel to hi. In contrast to the investi- 
gations of the Ising spin glass, none of these approaches 
have allowed to find numerically exact ground states with 
a reasonably high probability, such that, instead, effec- 
tively the properties of some set of metastable states with 
unclear relation to the ground-state behavior have been 
found and investigated. To improve on this, it is proposed 
here to combine exact ground-state computations of Ising 
variables embedded into the continuous spins with a spe- 
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cially tailored genetic algorithm exploiting the locally 
rigid cluster structure of metastable spin configurations 
|39l ] - This results in an efficient approach for ground- 
state computations of continuous-spin systems on planar 
graphs, which is tested and assessed here for the case of 
the bimodal XY spin glass on the square lattice, where 
numerically exact ground states can be found with high 
reliability for systems up to about 30 x 30 spins using cur- 
rently available computational resources. Implications of 
these results for the nature of the low-temperature phase 
of this model have been discussed elsewhere [H, Ho| ■ 

The rest of the paper is organized as follows. Section 
ITT1 is devoted to a description of the embedded match- 
ing technique for continuous spins, while in Sec. IIIII the 
combination of this approach with a genetic algorithm 
with cluster exchange is discussed. In Sec. IIV1 the per- 
formance of this approach for the 2D XY spin glass is 
investigated in terms of a detailed statistical analysis, fo- 
cusing on the large variations between disorder replica 
and offering a standardized approach of "quality assur- 
ance" for stochastic optimization algorithms. An exhaus- 
tive benchmarking of the new approach against general- 
purpose techniques is not feasible. At least, however, 
results comparing to the simple spin-quench used before 
and a more elaborate simulated annealing approach are 
presented in some detail. Finally, Sec. [V] contains my 
conclusions. 



II. EMBEDDED MATCHING 

For attacking the ground-state problem of continuous 
spin glasses in two dimensions, inspiration is taken from 
the polynomial solution of the Ising problem, which is 
hence described first, and then adapted to continuous 
spins via an embedding of Ising variables. 

A. Ising ground states as perfect matchings 

The polynomial complexity of the 2D Ising spin glass 
allows for the formulation of efficient algorithms for find- 
ing ground states and computing the partition function. 
A number of different techniques has been established for 
the calculation of the latter (nj 0, [H| , mostly relying 
on the computation of Pfaffians, but these will not be 
needed here. Computations of ground states rest on the 
concept of frustrated loops introduced by Toulouse [4l| : 
in the presence of couplings Jij of cither sign, for each 
closed curve along lattice links touching an odd number 
of antiferromagnetic bonds (Jij < 0), one cannot find 
a spin configuration satisfying all pair interactions, i.e., 
JijSiSj < for at least one ("broken") edge. Hence, 
the presence of such loops is responsible for the excess 
of the ground-state energy of the spin glass above the 
unfrustrated value -Efm = — Yliij) Due t° the con- 

tractibility of all loops on a planar graph, in this case it 
suffices to concentrate on the frustration of the plaque- 




FIG. 1: (Color online) Transformation of the Ising ground- 
state calculation on the square lattice to a matching prob- 
lem. Upper panel: Frustrated plaquettes (marked by small 
squares) have an odd number of antiferromagnetic (bold) 
bonds. The set of broken bonds forms a collection of lines on 
the dual lattice (shaded, gray lines). Lower panel: a ground 
state of the system is given by a minimum-weight perfect 
matching of frustrated plaquettes. The dashed line indicates 
an alternating cycle along which an exchange of matched and 
unmatched edges yields another perfect matching. 

ttes, i.e., the elementary faces of the lattice 12]. This 
is illustrated by the marking of frustrated plaquettes for 
the square lattice in Fig.rD For each plaquette □„, define 
the frustration function [4JJ 

$□„= Y[ signJy =±1, (3) 

such that <f>rj n = — 1 if and only if □„ is frustrated. By 
this definition, in a configuration of the Ising spins each 
frustrated plaquette must have an odd number (1 or 3 
for the square lattice) of broken bonds, whereas an un- 
frustrated plaquette is surrounded by an even number of 
broken bonds (0, 2 or 4 for the square lattice). Bonds 
drawn dual to the broken bonds then connect to form 
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energy strings starting and terminating in frustrated pla- 
quettes, cf. the sketch in the upper panel of Fig. [TJ The 
excess energy is 

i(£- J B FM ) = W s t ring = I J «I> ( 4 ) 

strings 

and a ground state corresponds to a collection of strings 
of minimum weight Wiring- Since, on a planar graph, 
closed loops of dual bonds can be contracted away, they 
cannot occur in a ground state, which hence corresponds 
to a minimum- weight perfect matching of the frustrated 
plaquettes. This is illustrated in the lower panel of Fig. [T] 
The planarity of the lattice ensures that each such match- 
ing corresponds to a valid spin configuration [1 21 ] . 

Following the above discussion, this matching problem 
is defined on the complete graph Q = T x T on the set 
T of frustrated plaquettes. Each of the \£\ = \T\\T — 1| 
edges e mn = (f m , /„) carries a weight 

VF(e mn )=min V" \J^\, (5) 

Tmn i ' 

(i,j)67m.™ 

corresponding to the minimum weight of all paths 7 mn 
on the (original) lattice connecting the plaquettes f m 
and /„. Hence an auxiliary minimum-cost path problem 
needs to be solved as an input to the matching calcula- 
tion. This is most efficiently achieved by an appropriate 
implementation of Dijkstra's algorithm with 0(|£| In \£\) 
complexity, or, for the case of a bimodal P(Jjj) where 
\Ja \ = Jo for all edges, by a simple breadth-first search 
[421 ] . Since there is an even number of frustrated plaque- 
ttes [78{, a perfect matching on Q can always be found. A 
polynomial algorithm for the matching problem on gen- 
eral graphs has been proposed by Edmonds [431 ] . It pro- 
ceeds by successively identifying augmenting paths in the 
matching graph, i.e., cycles of alternating matched and 
unmatched edges such that an exchange matched <-> un- 
matched decreases the overall weight. This is illustrated 
by a cycle in the original lattice in the lower panel of 
Fig. [U The complexity of the original implementation 
is O ( | 2 1 ^ | ) . The present implementation is based on 
the "Blossom IV" matching algorithm of Cook and Rohe 
incorporating many improvements developed in the com- 
binatorial optimization literature after Edmonds' original 
proposal [441 ] . 

Given a solution to the matching problem, a corre- 
sponding spin configuration is found by arbitrarily choos- 
ing the orientation of one spin and successively imple- 
menting the satisfaction constraints expressed by the per- 
fect matching via selecting spin orientations in a breadth- 
first search emanating from the chosen starting point. 
Note that, depending on the distribution of couplings 
Jij , neither the solution of the matching problem nor the 
map ping back to spin configurations needs to be unique 
[4a |: if some edges in the matching problem have the 
same weight, there could be different matchings of min- 
imum weight. On the other hand, also the notion of 
minimum-weight paths on the original square lattice is 




FIG. 2: (Color online) Embedding of Ising spins into the con- 
tinuous rotators Si via decomposition with respect to a di- 
rection r in spin space. 

not necessarily unique, and more than one path between 
two frustrated plaquettes could be of minimal weight. Fi- 
nally, each configuration of energy strings corresponds to 
two different spin states, related by spin inversion. This 
generally leads to a large ground-state degeneracy for dis- 
crete and rational distributions P{{Jij}), but a unique 
ground state, e.g., for the Gaussian case [46[. For the 
complete graph Q, the time complexity of Edmonds' im- 
plementation would be 0(|J-| 4 ), corresponding to 0(L 8 ) 
for a L x L lattice. Although this is polynomial, further 
improvements are highly desirable to reduce the rather 
large exponent and enable treatment of reasonably sized 
problem instances. This can be achieved by a thinning of 
the complete graph Q: a matching of frustrated plaque- 
ttes in the ground state becomes more and more unlikely 
with increasing weight of the path connecting them, and 
such large- weight edges can consequently be disregarded. 
Suitable cutoff parameters depend on the distribution 
of Jy, and have to be tested thoroughly. For the im- 
plementation used here, cutoffs at fixed maximum path 
weight and conditions on the minimum vertex degree in 
the matching graph have been employed with comparable 
success. 



B. Embedded matching for continuous spins 

For continuous spins, the notion of plaquette frustra- 
tion stays meaningful, since it is a property of the bond 
configuration only. The transformation to a matching 
problem, however, is restricted to discrete Ising spins. To 
leverage the tractability of the Ising case for the treat- 
ment of continuous-spin systems, an embedding of Ising 
spins into the continuous rotators is employed. To this 
end, consider an arbitrary direction r, |r| = I, in spin 
space common to all lattice sites, and decompose the 
O(n) spins Si of Eq. Q as Si = S l }+S± = (S z -r)r + S^, 
cf. the illustration in Fig. [2] This induces a decomposi- 
tion H = H r > 11 + U r ^ with 

H r,\\ = _J2jr. e r e r p (6) 
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where = sign(S , . i 
J[j are given by 



r), and where the effective couplings 



0.3 



J ij | 



(7) 



Hence, with respect to reflections of the spins Si along 
the plane defined by r, the signs e[ are Ising variables 
(cf. Fig. [2]), and the embedded Hamiltonian ((6| is that of 
an Ising model. Note that with respect to these reflec- 
tions Si i— > Si — 2(Si-r)r, the perpendicular part H r ' ± is 
invariant and thus does not contribute to the embedded 
dynamics. A similar embedding of Ising variables has 
been used to formulate a cluster-update Monte Carlo al- 
gorithm for continuous spins |47| . 

Updating the effective Ising variables e£ i— > — e\ via 
plane reflections Si i— > Si — 2(S< • r)r, a ground state of 
the Hamiltonian 7i r '" of ([6|) can be found, for instance, 
using the transformation to a matching problem outlined 
above. Note that since sign J^- = sign Jy according to 
Eq. (|7|), the frustration function <&□ does not depend 
on the embedding direction r, and only the weights of 
the matching graph Q must be updated for each embed- 
ded matching computation. The rotational symmetry of 
the O(n) Hamiltonian |T]) can then be recovered by a 
random sampling over different embedding directions r. 
This leads to the following algorithm: 



1 
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procedure EmbeddedMatching({S' 1 }, J) 
for j <— 1,1 do 

choose random direction re- 
determine ground state of W J '"(e[ J ) 
for i <— 1, L 2 do 
if e^ 3 ^ e[ J then 

Si < Si — 2(Si ■ fj)fj 
end if 
end for 
end for 
end procedure 



For each direction r, it holds that H = W'"({el}) + 
H r.± > H r 'H({e[}) +U r ' X , where e\ is the ground state 
configuration of the embedded Ising model. Conse- 
quently, the embedded matching procedure corresponds 
to a strictly downhill minimization approach. If TL is in a 
ground state, Ti r '" ({e[ }) = TL — TL 7%± must be in a ground 
state for each r as well. Conversely, however, H r '"({e[}) 
being in a ground state for each r does not guarantee 
global minimum energy for the full TL. This is due to the 
fact that the embedded couplings of depend on 
the spin configuration {Si} and hence on the history of 
previous embedding matching runs. As a consequence, 
the dynamics of {Si} induced by the embedded match- 
ing procedure has metastable states [7!| . The number of 
metastable states is far less, however, than for the local 
spin-quench approach of Eq. (J5J), since for each direction 
r a global minimum is found: while TL Ti ^ converges to a 
locally spin-flip stable state (otherwise a direction r could 
be found, for which embedded Ising minimization would 
lead to reflections of one or more spins), not every such 
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FIG. 3: (Color online) Histograms of energies of metastable 
states of a 24 x 24 sample of the ±J XY spin glass on the 
square lattice as found from 1 000 independent runs of the 
embedded matching technique compared to results of a lo- 
cal spin quench. The onset of the hatched area to the right 
indicates the ground-state energy of this sample. 



locally stable state is metastable with respect to the em- 
bedded matching procedure (because the embedded Ising 
system for some direction r could be metastable instead 
of globally optimal). 

It is found numerically that the sequence {Ei} of en- 
ergies of the embedded matching approach always con- 
verges. As a consequence of the dependency of of 
([7]) on {S^}, the limit depends on the particular se- 
quence {ri, r 2 , . . .} of chosen directions. Figure[3]shows a 
histogram of energies found from the embedded matching 
approach with a number of different random starting con- 
figurations and different series of embedding directions 
for a particular 24 x 24 sample of the 2D bimodal XY spin 
glass. For comparison, a corresponding histogram for the 
local spin-quench method of Eq. (|2|) is also shown. It is 
apparent that the average energy of metastable states 
is lower for the embedded matching technique, and this 
behavior is found to survive averaging over the random 
couplings {Jij}- Nevertheless, the probability for con- 
verging to a ground state is apparently very small for the 
system size considered, cf. Fig. [3] One might speculate 
that this shortcoming is connected to the fact that only 
reflections of spins along a plane, i.e., improper rotations, 
are allowed updates in this approach: if the intermedi- 
ate, improperly rotated configurations connecting a state 
to another, properly rotated state of lower energy have 
higher energy, they form a barrier which cannot be over- 
come by a strictly downhill procedure. The only other 
transformation 1Z besides plane reflections allowing for 
an Ising type symmetry 1Z 2 = id are point inflections 
Si i—* —Si. The embedded matching technique can be 
extended to include these transformations. Their inclu- 
sion, however, is not found to remove significantly many 
barriers, such that this approach is not further considered 
here M. 
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III. BOND-ENERGY DIFFERENCE 
CROSSOVER AND GENETIC MATCHING 

Although an important improvement over the local 
spin quench approach ^ty, embedded matching alone is 
not sufficient for reliably finding ground states. Further 
advances are possible by an understanding of the struc- 
ture of metastable states exploited in a suitably tailored 
genetic algorithm. 



A. Rigidity and domain structure 

To understand the mechanism of metastability in the 
embedded matching approach and develop a strategy for 
overcoming it, one needs to take into account some fea- 
tures of the low-temperature phase of spin glasses. While 
there is no long-range order, the freezing of spin orien- 
tations corresponds to some short-range order, expressed 
in a non-zero range of correlations • Consequently, at 
low temperatures spins are rather rigidly locked together 
locally, and their orientation can only be changed (at 
low, but generally non-zero energies) by a rigid O(n) ro- 
tation of a cluster of spins J48|. Therefore, the manifold 
of internal states (i.e., the parameter space of the rele- 
vant order parameter) is described by the full orthogonal 
group O(n), in contrast to the case of homogeneous mag- 
nets, where the global magnetization confines the internal 
states to the quotient space SO(n)/SO(n — 1) ~ S n , i.e., 
an n-dimensional unit sphere [iij IHJ. Such spin clus- 
ters hence behave like solid n-dimensional bodies. Note, 
however, that their O(n) rotation is not in general a 
zero mode, but a low-energy excitation. The existence 
of such clusters could recently be explicitly revealed uti- 
lizing the genetic embedded matchi ng appr oach for the 
planar spin glass in two dimensions [33 . |4Q | . This sym- 
metry also determines the topologically stable defects in 
spin glasses: as in ferromagnets, they are determined by 
the homotopy groups of the internal space, here 0(n). 
For planar rotators, for instance, in addition to vor- 
tices (resp. vortex lines) also present in the homogeneous 
case, this framework predicts domain walls, which can 
be directly observed in form of chiral walls for the (two- 
dimensional) XY spin glass poL[5lj|. Consequently, some 
important classes of low-energy excitations in continuous 
spin glasses are: 

1. Rigid O(n) rotations of spin domains. 

2. Topological defects: domain walls, vortices etc. 

3. Smooth, spin-wave excitations. 

In the context of ground-state searches, spin waves can 
be easily removed by local relaxation techniques (see the 
discussion below in Sec IIIIB| ). Some of the topological 
defects, such as domain walls, can be composed out of a 
sequence of domain rotations, such that I concentrate on 
this first type of excitation here. Note that the given clas- 
sification is not meant to be exhaustive, i.e., it does not 
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FIG. 4: (Color online) Relative domain decomposition of 
metastable states of the two-dimensional XY spin glass, (a) 
Density plot of bond energy differences ([8j for two metastable 
states of the embedded matching method for a given 24 x 24 
disorder sample, (b) Cluster decomposition resulting from 
the Hoshen-Kopelman algorithm with clustering rule (O and 



cutoff parameter k 



0.3a 



a,/3 



0.05. 



53], 



express a prejudice as to whether the asymptotic low- 
energy excitations in spin glasses are droplets [52 
mean-field like extended defects [Hi[ or "sponges" 
for instance. It is, instead, only used as a guideline for 
the identification of appropriate metavariables in the for- 
mulation of an efficient ground-state search heuristic. 

Given that typical metastable states differ by rigid 
O(n) rotations of domains, an explicit implementation 
of such rotations in an optimization heuristic enables it 
to perform a search directly on the space of metastable 
states. Direct inspection of the transformations con- 
necting metastable states in continuous-spin glasses show 
that this domain structure indeed is a valid description, 
see Refs. |4Cj, |5l|, |56(. Note that the concept of such do- 
mains is a relative one, i.e., the domain decomposition 
of a metastable configuration can only be determined 
with respect to other metastable states. In particular, 
a domain decomposition may be performed for a pair of 
configurations, identifying the set of O(n) rotations map- 
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ping them onto each other. This might be achieved by 
determining, by a singular-value decomposition, locally 
averaged rotation matrices connecting the configurations 
[40L l5a]. Here, a (computationally) simpler approach is 
chosen by noting that the O(n) symmetry of the Hamil- 
tonian ([1]) ensures that for each such domain all energies 



J ij Si 



Sj of bonds in the interior are invariant, 



whereas the energies carried by bonds crossing the do- 
main boundary change due to local mismatches of the 
surface spins of the rotated domain and its environment. 
Some of these changes will be of spin- wave type and hence 
will have been relaxed away once metastability has been 
reached again. Some differences, however, remain, giving 
a handle on domain identification. Consequently, for a 
pair (a, f!) of metastable configurations one might con- 
sider the bond energy differences (BEDs) 



Jij {Sf ■ Sf 



(8) 



The distribution of BEDs is depicted for a pair of 
metastable states of the 2D XY spin glass in Fig. SJa), 
showing clear structures of rigid domains. Defining a do- 
main boundary by a BED exceeding a threshold value, 



i.e., 



\ae; 



> K 



(9) 



a domain decomposition can be performed, for instance, 
with the Hoshen-Kopelman algorithm [57|. This is il- 
lustrated in Fig. IHb) for the BEDs of Fig. The 
cutoff parameter K a '@ is chosen here of the order of the 
total variation of BEDs over the whole sample, i.e., pro- 
portional to the standard deviation cr^ D , to accommo- 
date differences between disorder realizations as well as 
metastable states of largely varying energies. Since these 
domains are merely utilized for a more efficient ground- 
state search, I do not have to bother here with the ques- 
tion of whether there is a precisely defined characteristic 
length associated with such domains [56} , independent of 
the size of the system, as having k 



a,/3 



a,0 
r BED 



detects 



the length(s) appropriate to the sample at hand auto- 
matically. 



B. Genetic matching 

Embedded matching in combination with domain de- 
composition of configurations with BED clustering al- 
lows optimization directly on the space of metastable 
states of the embedded matching method. This already 
corresponds to an enormous reduction in the size of 
the phase space. The meta-optimization on metastable 
states is performed here utilizing a hybrid genetic algo- 
rithm, although variants based on other global optimiza- 
tion strategies such as simulated annealing are conceiv- 
able as well. Genetic algorithms [581 ] mimic natural evolu- 
tion by maintaining a population of candidate solutions, 
which is evolved in generations by a process involving the 



crossover and mutation of solutions followed by a selec- 
tion of members with higher fitness, i.e., lower energy for 
the case of ground-state computations considered here. 
In the canonical form of genetic algorithm, solutions are 
represented by bit strings in a binary representation and 
crossover and mutation correspond to the random ex- 
change of bits between solutions and random bit flips, 
respectively (58j . In this form, genetic algorithms have 
been applied to the Ising spin glass, but, unless for very 
small systems, true grou nd states could not be found 
with high reliability [59l [60j . Only hybrids combining 
genetic crossover with some downhill optimization proce- 
dure such as local spin flips or the "cluster-exact approx- 
imation" [l9j], restricting the search space to metastable 
states, led to more successful approaches [20I [olj. 

Although widely and successfully employed, there is 
no theoretically sound framework for designing efficient 
genetic algorithms [58| , such that their construction rests 
on heuristic strategies and additional insight specific to 
the problem at hand. Generally, one strives to achieve 
a balance between fast convergence to an optimum an- 
swer and the upholding of genetic diversity throughout 
the "evolution" (which, in turn, tends to slow down con- 
vergence), in order not to miss the global optimum. The 
choice of crossover operation appears to be most impor- 
tant in this context. In the present work, instead of ran- 
domly exchanging single spins, the BED cluster decom- 
position is employed to exchange domains of rigid spins 
between solutions. This has the important advantage of 
retaining the high degree of optimization already found 
from the embedded matching technique inside of the do- 
mains and directly operating on the space of variables 
relevant for the construction of metastable states. The 
domain decomposition can be performed directly with 
the "parent" configurations to be combined ("diadic" 
crossover) or, alternatively, by using a third, "mask" con- 
figuration from the population used only for the domain 
decomposition. The latter, "triadic" crossover is used 
here, similar to the technique suggested in Ref. [62j for 
Ising spins, since it is found to perform slightly better 
for continuous spins. Genetic diversity is strengthened 
by restricting the selection of parents to be combined to 
neighbors after the population has been arranged in a lin- 
ear ring [62j]. This introduces some degree of geometric 
"locality" of the population and allows good solutions to 
be refined independently in different areas of the config- 
uration set. Previous approaches [2(| [fill [H[ have used 
a fixed total number of crossover operations per mem- 
ber of the initial population, followed by a halving of the 
population by elimination of the higher-energy solution 
of each pair of neighboring configurations, and then a re- 
iteration of the remaining population. This reduces the 
total effort by removing unpromising solutions and bring- 
ing distant parts of the "ring" of solutions closer to each 
other in later stages of the optimization. For the time 
being, I will adopt this technique here as well. A more 
efficient variant, geared at the detection of hard samples, 
is presented below in Sec. IIV CI In total, the resulting 
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genetic embedded matching (GEM) algorithm proceeds 
as follows: 



notes the corresponding global rotation matrix, and 



a/3 _ J_ \ " na a/3 
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procedure GEM(S, C, V,I, C) 

initialize { } , k — 1 , . . . , S randomly 
s <- «S 

while s > 4 do 

for c <— 1 , C x s do 

randomly select pair (a, /3 = a+1) of confs. 
randomly select mask conf. 7 ^ a, 
{Sf,Sf} = BEDCROSSOVER(a, /3, 7) 

MUTATE({Sf},7>), MUTATE({Sf },7>) 
EMBEDDEDM ATCHING ( { Sf } , I) 
EMBEDDEDMATCHING({Sf 

for j <— 1 , £ do 
for i <— 1, L 2 do 

S? «- hf /|hf 1 

end for 
end for 

if }) < H{{S? }) then 

S? ^ S?,i = l,...,L 2 
end if 

ifW({Sf})<W({Sf}) then 

Sf^§?,i = l,...,L* 
end if 
end for 

for all distinct pairs (a, a + 1) do 
if H({S?})<H({S?+ 1 }) then 

remove configuration a + 1 
else 

remove configuration a 
end if 
s <— s — 1 
end for 
end while 

output (best of) remaining configurations 
end procedure 



(10) 



is the matrix of overlaps. This maximization is performed 
by a singular value decomposition to diagonalize , in 
which case q a ^ is just the trace of the resulting diagonal 
Qa? 0, HH ■ This form of replacement restriction helps 
to maximize genetic diversity [2(| ■ After sC crossovers, 
the higher energy instance of each adjacent pair of config- 
urations is discarded, thus halving the population. The 
complete process is repeated until at most four configu- 
rations arc left, which form the result of a run. 

As will be shown in the next Section, this combina- 
tion of techniques in the genetic embedded matching 
method allows for the determination of (numerically) ex- 
act ground-states of reasonably large continuous-spin sys- 
tems in 2D with high reliability. 



IV. PERFORMANCE 

Using probabilistic methods for ground-state searches, 
special care is needed to ensure that true ground states 
are found. Since for MV hard optimization problems 
the decision variant is MV complete, there is no way 
of definitely distinguishing a metastable from a ground 
state short of an exact solution of the instance. A gen- 
eral probabilistic approach of "quality assurance" for the 
GEM method is outlined and applied to the 2D XY spin 
glass in Sec. IIVBI In some dynamical approaches, such 
as local spin- flip Monte Carlo simulations, the specific 
hardness of a sample shows up in the behavior of auto- 
correlation times, to which a simulation run can in prin- 
ciple react dynamically by increasing the simulation time 
accordingly. Below in Sec. IIVC1 it is discussed whether 
a similar heuristic for detecting hard samples can be ap- 
plied in the GEM approach. 



The BEDCrossover operation performs the BED 
domain decomposition of Sec. IIII Al with respect to the 
mask and, for each domain, either swaps all spins be- 
tween the parents a and j3, copies domains a 1— ► (3 or 
(3 1 y a, or leaves the domain invariant, with all possi- 
bilities occurring randomly at equal proportions. Muta- 
tions are performed by randomly choosing new spin ori- 
entations with a probability V. The resulting offspring 
are optimized using J iterations of EmbeddedMatch- 
ING from Sec. Ill BL followed by C iterations of a local 
spin quench. The latter is useful since, close to a mini- 
mum where only spin- wave excitations are left, both ap- 
proaches converge to the same state, but the spin quench 
is much faster computationally. Lower-energy offspring 
then replace their parents. In the implementation used, 
each offspring is only compared to the morphologically 
closer parent, i.e., the one with a larger optimized scalar 
overlap q al3 = max^o^) E ff , T 1%? R %? > where R *? de ~ 



A. Performance and comparison to simulated 
annealing 

Local spin quenches according to Eq. ([2| yield states 
in a broad range of energies, cf. Fig. [3] For ascribing the 
ability to find ground states to a stochastic method one 
would, instead, require that states of exactly the same 
energy (up to machine precision) are found in a sizable 
fraction p s of attempts (withp s = 95%, for instance) and 
that no states of lower energy can be found with runs 
of largely increased effort or utilizing other optimization 
techniques. As is evident from Fig.[3l this also cannot be 
said of the embedded matching approach alone. Figure[S] 
shows the minimum energies found in repeated runs of 
the GEM technique for the bimodal XY spin glass in 
2D with a randomly picked disorder realization of size 
24x24 (which is identical to the realization used in Fig. [3]) 
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FIG. 5: (Color online) Histograms of energies found for a 
24 x 24 sample of the bimodal XY spin glass from 100 runs of 
the genetic embedded matching approach with a population 
of size S = 64 (hatched bars, average runtime 700 s on a 
Pentium IV 2.8 GHz) as compared to simulated annealing 
runs with a total of about 2 x 10 7 Monte Carlo sweeps (solid 
bars, average runtime 4 000 s). The disorder realization is 
identical to the one considered in Fig. [3] 



and a starting population size S — 64. For comparison, 
this Figure also shows the histogram of repeated runs 
of an extensive simulated annealing [H, [gl] computation 
with an exponential temperature protocol (a linear pro- 
tocol yields very similar results) and a total number of 
about 2 x 10 7 lattice sweeps of single spin flips per run, 
leading to an about sixfold runtime as compared to the 
GEM computations. As is seen, the GEM runs result 
in clearly lower energies than the simulated annealing. 
Additionally, the latter still show a sizable spread of the 
energies found, whereas the GEM technique appears to 
yield states of the same energy. Only on going to much 
higher energy resolution, the GEM results are resolved 
into a small number of distinct peaks cumulated from 
runs yielding identical energy up to (or close to) machine 
precision, cf. the inset of Fig. \5\ A fourfold increase of 
the starting population to S — 256 leads to a convergence 
of all 100 runs to the lowest-energy peak to the right in 
the inset of Fig. [5] No further increases of the popula- 
tion size up to S = 1024 lead to lower energies such that, 
with high confidence, this peak corresponds to the true 
ground-state energy of the system. Consequently, it can 
be said that runs with S — 64 have a probability of about 
p s = 16% of leading to a ground-state. Methods for guar- 
anteeing high reliability of finding ground states over the 
distribution of disorder realizations are discussed below 
in Sees. HVBl and ITTCl 

Although, for the given disorder realization, the GEM 
technique appears able to find ground states and clearly 
outperforms the simulated annealing approach, varia- 
tions in the "hardness" of different replica in the random 
couplings are known to be large (see, e.g., Refs. [65l.[66T|). 
and the corresponding variation in the efficiency of the 
methods should be taken into account. Since the com- 



FIG. 6: (Color online) Average minimum energy {E(T))j for 
1 000 samples of size 16 x 16 of the bimodal 2D XY spin glass 
found from GEM and simulated annealing runs with a total 
runtime T (in re-scaled units). 



monly considered distributions P(Jij) contain the ferro- 
magnetic lattice with Jy = Jo > 0, which is trivially 
handled by either optimization method, the behavior of 
interest can only be either that for the worst case, which 
is, however, difficult to assess for the spin-glass model 
considered, or rather the average performance for the dis- 
order distribution at hand. As a first step in this analysis, 
I considered the convergence of the average minimum en- 
ergy observed with the computational effort invested. For 
simulated annealing with Metropolis acceptance rule, it 
is known that with logarithmically slow cooling, ground 
states will be found in finite (but, of course, very large) 
time 
ever. 



67[. Since this cooling schedule is impractical, how- 



exponential or power-law cooling curves are used 
instead [HI, |64| ■ The possibility of different acceptance 
rules complicates things further, and it is naturally im- 
possible to benchmark against all these variants. I re- 
strict myself here to the probably most commonly used 
exponential protocol. The asymptotic form of energy 
convergence in simulated annealing of spin-glass systems 
has been the topic of some debate in the past 0, [H|. 
Numerically, a power-law convergence, 



A F T 



-C 



11) 



for large cooling times T (i.e., the total number of Monte 
Carlo sweeps) was found for the 2D Ising spin glass, 
while, on the contrary, logarithmic convergence was ob- 
served for the 3D variant [25| ■ On the other hand, on the 
basis of modeling (Ising) spin glasses as sets of weakly in- 
teracting two-level systems, it was conjectured that the 
logarithmic form would be universal [5J. Here, (-)j de- 
notes the average over disorder. Figure [5] shows the aver- 
age minimum energy found from simulated annealing of 
1 000 samples of size 16 x 16 of the 2D XY spin glass for 
annealing times between T = 50 000 and T = 3.2 x 10 6 
sweeps, compared to the energies found from GEM runs 
with population sizes S oc T between 8 and 128 configura- 
tions. The abscissa for the simulated annealing data has 
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been linearly rescaled to result in equal runtimes for both 
approaches on a Pentium IV 2.8 GHz (both algorithms 
scale linearly in T). The data from simulated anneal- 
ing can be fitted with reasonable quality to the power- 
law form (fTTj) . yielding a decay exponent £ = 0.33(36), 
whereas a logarithmic form does not adequately describe 
the data. This is comparable to the value C = 0.25 
found for the 2D Ising spin glass in Ref. [25|]. Note 
that the extrapolated asymptotic ground-state energy 
Eva = —387.45(76) is compatible statistically with the 
value found from GEM runs already for the smallest pop- 
ulation size 5 = 8 considered. In fact, the GEM data are 
constant within statistical errors for S > 32. The varia- 
tion of energies found from the GEM technique can also 
be described by (fTTj) . resulting in £ ~ 2.3(48) (the large 
statistical error results from the only minute variation of 
(E) observed). 

The GEM algorithm as presented in Sec. IIIIBI in- 
volves a number of parameters which need to be tuned 
to achieve these good results. Performance appears to be 
rather weakly dependent on the mutation rate, and best 
results are found for a rate of about V — 2.5%. Much 
more frequent mutation destroys the relatively good op- 
timization achieved at intermediate stages and decrease 
the overall performance. Since the offspring configura- 
tions produced by the BED crossover are still optimal 
inside of domains, relatively small numbers of embedded 
matching and local relaxation steps are found to be suf- 
ficient, T — 15 and C = 100 was usually chosen here 
(cf. the pseudocode of the algorithm in Sec. IIIIB]) . The 
number C of crossovers per replica determines how well 
the available "genetic pool" of original configurations is 
explored. Beyond a certain number of crossovers, the 
population becomes uniform and further increases do not 
improve the probability of finding ground states. For ac- 
cessible system sizes C = 8-16 is a good choice. The main 
tuning parameter of the approach is the initial popula- 
tion size S, which is changed to accommodate the vari- 
able hardness of different system sizes, models and disor- 
der realizations. It is the only parameter whose increase 
ultimately guarantees ground states to be found. Note 
that the total number of crossovers is 2C{S— 4) (assuming 
S = 2") and hence linear in S. For a given single realiza- 
tion, computation of true ground states can be guaran- 
teed with high confidence by tackling the same disorder 
configuration with largely increased computational effort 
(in particular, say, a fourfold increase in population size 
S), until no further change in energy is observed. For the 
random distribution of configurations to be investigated, 
however, a more automatic (and less computationally ex- 
pensive) approach is required. 



B. Probabilistics of successes 

For stochastic optimization methods, arrival at true 
ground states cannot be guaranteed. For given input 
data in form of the disorder realization and a choice for 




FIG. 7: (Color online) Estimated failure probability pf — 1 — 
p a for the GEM technique applied to XV spin-glass samples 
of size 16 x 16 as a function of the computational effort T — S. 
The lines show fits of the form (|12p (easy and hard sample) 
resp. (TT3J (disorder average) to the data. 



the tunable parameters, a ground state is found with the 
success probability p s ({ Jy}; T), where T denotes the rel- 
evant parameters. As was discussed above, by far the 
most influential parameter for the GEM approach is the 
initial population size S, such that I here restrict consid- 
erations to T = {S}. Full information about the distri- 
bution of p s induced by P{{Jij}) and the dependency on 
the algorithm's parameters would correspond to a com- 
plete understanding of the performance characteristics or 
generalized computational complexity Q • This computa- 
tion, however, is impractical due to the high-dimensional 
nature of this parameter space: using, e.g., 100 runs to 
estimate p s for a given set of parameters for 1 000 disorder 
realizations and 100 combinations of parameters S,I, . . . 
would require 10 7 ground-state computations for a sin- 
gle system size! From p s ({ </ij}; T) one could deduce the 
perhaps most interesting distribution T m j n ({ Jy};p s ) of 
efforts required for a constant success probability p s . 

Figure [7] shows the failure probabilities pf = 1 — p s for 
an "easy" and a "hard" sample of the 2D bimodal XY 
spin glass as a function of the population size S = T, 
compared to the average failure rate pf over 100 dis- 
order replica. The huge spread in hardness is apparent: 
while, for instance, only in 2 out of 100 cases do runs with 
S = 64 fail to find a ground state for the easy sample, 
74% of attempts for the hard sample end in a metastable 
state. It is therefore not enough to fix the run parame- 
ters by considering one or two randomly chosen configu- 
rations. For a description of the functional form of Pf(T) 
in Fig. [71 consider performing n statistically independent 
runs of length To with failure probability pf$ and picking 
the solution of lowest energy as final answer. With this 
prescription, a ground state is not being found only if all 
of the runs fail, and the combined failure probability is 
thus 

Pf (T = nT )= P T f / To . (12) 



11 



lO" 2 f=r 



T = S = 


12 






T = S = 


64 




1(T 3 


-T = S = 


192 






g io- 4 

10" 



10" 




500 1000 1500 2000 



FIG. 8: (Color online) Histogram H(pf) of failure probabili- 
ties for 1 000 disorder samples for the 2D XY spin glass as a 
function of initial population size S. The lines show analytical 
approximations discussed below in the main text. 



FIG. 9: (Color online) Estimated probability-density function 
G(2mi n ) of minimum required runtimes T m i n for 16 x 16 sys- 
tems of the 2D XY spin glass to achieve success probability 
p s = 95%. The dashed line shows a fit of the generalized 
extreme-value distribution (|15p to the data. 



Due to the locality constraint in choosing parent config- 
urations for crossover, increasing the initial population 
size S = T has essentially the same effect as performing 
independent runs. For sufficiently large T, Eq. (fT^J) is 
hence an excellent description of Pf(T) for a single sam- 
ple. Figure [7] shows fits of the form (TPS]) to the data for 
the "easy" and "hard" samples. There is only a single fit 
parameter, p 1 /^ , corresponding to a measure of hardness 
of the sample (with respect to the GEM technique). The 
ensemble average {pf({Jij}', T)} j cannot be expected to 
follow the same exponential form (fl"2"|) since, in general, 
(exp[aT]} ^ exp[(a)T]. It is found, however, that it is 
well described by the slight generalization 



(P/({J«};T)). 



A n T/To 



(13) 



with an additional amplitude A p < 1, as is apparent 
from the corresponding fit also shown in Fig. [7] Con- 
sequently, the average failure probability decreases more 
slowly with T than would be expected from the behav- 
ior on single samples. Note that due to the form l|12p 
it is not appropriate to consider the combination T/p s 
as the "computational effort" of a sample [66| , since this 
assumes a linear relation between p s and T. 

In view of the results of Fig. [7J it is of interest to inves- 
tigate the distribution H(j>f) of failure probabilities over 
disorder samples. To this end, the failure probability was 
sampled by performing 100 independent runs for each of 
1 000 disorder samples with run length T = S = 64. 
The histogram estimating the probability-density func- 
tion H(pf) is shown in Fig. [8l revealing clearly the 
breadth of this distribution, reflecting the large spread in 
hardness already suggested by the data of Fig. [7J From 
the histogram H(j>f), it is possible by means of Eq. (TT^]) 
to recover the distribution of the minimum required run- 
times T m i n ({ Jij};pf), corresponding to the distribution 
of hardness of samples under the GEM technique: from 
the estimate of Pf({Jij}) for a disorder sample {Jij} at 



fixed runtime T, Eq. (fT2]) implies 



r min ({j lj }) = T 



lnp f)( 



]n.p f ({Jij}) 



(14) 



Figure [9] shows the distribution G(T m i n ) of minimum 
runtimes at failure probability ptfi = 0.05 thus result- 
ing from the runs at fixed runtime T = 64 presented 
in Fig. [8l The breadth of the distribution is apparent: 
while the average is around (T m i n }j ~ 225, there is a fat 
tail with some of the 1 000 disorder configurations fea- 
turing a Tmin as large as 2 000. Such density functions 
typically occur when considering the extrema of large 
samples drawn from underlying probability distributions. 
Asymptotically, the distribution of extremes is known to 
be universal, following the form [68[ 



x exp 



x — fl 



-1/?' 



(15) 



where the parameter £ depends on the tail behavior 
of the underlying, primary distribution for large argu- 
ments x. Depending on £, this form is known as Weibull 
(£ < 0), Gumbel (£ -> 0) or Frechet (£ > 0) distri- 
bution, respectively. As is seen in Fig. [9l it fits the 
data for T m i n extremely well, resulting in £ = 0.270(44), 
H = 115.6(45), and a = 87.7(39) with an excellent 
quality-of-fit Q = 0.23. Similar distributions of the 
Frechet type have been found for the tunneling times in 
Wang-Landau flat-histogram simulations of the Ising spin 
glass [66J, |69j . One might speculate on the origin of the 
occurrence of extreme-value statistics in hardness mea- 
sures of spin-glass samples: if, as has been suggested [f|, 
a spin-glass sample can be described as a set of n — n(L) 
weakly interacting two-level systems, it appears plausi- 
ble that the largest barrier or the slowest relaxation time 
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determines the hardness of the configuration. Then, the 
hardness would be the maximum or minimum of a sam- 
ple of size n from the underlying distribution of two- 
level systems, asymptotically distributed according to the 
extreme- value distribution (| 1 5[) . In line with this argu- 
ment, it was recently observed [70( that the distribution 
of relevant barriers in the Sherrington-Kirkpatrick model 
follows a Frechet distribution with a value of £ ~ 0.33, 
rather similar to the form found here. 

Given that G^ i/JiCr (T m i n ) describes the GEM data so 
well, it is worthwhile to use Eq. (| 1 2|) to reveal the result- 
ing analytical form of the distribution of failure probabil- 
ities, H(pj). Following standard probability theory fn\ . 
density functions transform as 



H(p f )dp f = G^ a (T u 
which, using Eq. (fT2|) . leads to 



dT, T 



dp f 



■dp/, 



p f (In p f ) 2 In p f 



(16) 



(17) 



As is seen from the curves in Fig. [8j this form with the 
parameter values £, /i and a given above fits the numer- 
ical distribution H(pf) for the same T = 64 data per- 
fectly well (i.e., the approach is self-consistent). Addi- 
tionally, however, it describes independent runs of differ- 
ent lengths to high precision and hence the form (|17p is an 
excellent description for general runtimes T, as indicated 
by the additional curves in Fig. [5] Consequently, the 
three-parameter family of distributions (| 1 5|) and the lim- 
iting distributions derived via Eq. (fT2"|) form a complete 
description of the full probability density p/({Jij};T). 

While it certainly would be instructive to extend the 
analysis of T min via the distributions (fl~5)) to a scaling 
analysis of the fit parameters £, fi and a with lattice 
size L, the huge computational effort would be dispro- 
portional. Instead, I concentrate on the mean required 
effort (T m i n ) j as a function of system size, averaged over a 
smaller disorder sample of only 100 configurations. These 
data for failure probability pf — 5% are shown in Fig.fTUl 
together with a fit to the expected exponential form 



(T m in) j — A T e 



(L/L ) 



(18) 



which works reasonably well with parameters At — 
4.82(26) and L = 8.515(89). Increasing the rate of toler- 
ated failures to, e.g., Pf = 10%, merely reduces the pre- 
factor to At — 3.64(19), but leaves Lq almost invariant. 
This data brings back to attention the fact that, although 
the GEM approach works quite well, and clearly outper- 
forms simulated annealing, it naturally cannot evade the 
AA'P-hard nature of the problem enforcing an exponential 
growth of computational effort. To complicate the matter 
further, it is well conceivable that the shape parameter 
£ of (|15|) increases with system size, as was observed in 
tunneling simulations of spin-glass models [66] . Since for 
£ > 1/2 the variance of G^ iAliCr becomes ill-defined, and 




FIG. 10: (Color online) Average effort (T m i n )j for finding 
ground states with constant success probability p 3 = 95% 
as a function of lattice size L 2 , estimated from 100 disorder 
configurations per size. The line shows a fit of the form (|18[) 
to the data. 



for £ > 1 additionally the mean diverges, this would im- 
ply that a correct choice of population size S = T for 
all disorder configurations becomes impossible beyond a 
certain system size. 



C. Hardness of samples and refinements 

Numerical ground-state computations in (non- 
polynomial) spin-glass systems are subject to two types 
of "hardness problems" : the exponential growth of aver- 
age computational effort with system size and the large 
fluctuations in sample hardness implied by heavy-tailed 
distributions of the type (|15p . While AA'P-hardness 
means exponential effort for the worst-case samples, it 
is clear that (close to) ferromagnetic (i.e., best case) 
configurations can be tackled in polynomial time. Hence, 
the difference in effort diverges with system size. While 
this is true for the set of all possible input data, it is not 
clear a priori that the samples receiving non-negligible 
weight from the P({Jij}) considered indeed show such 
spread as implied by the data of Fig. In this Section, 
I discuss technical refinements of the GEM technique 
designed to address the problem of large fluctuations in 
hardness. 



1. Effort adaptation 

The fixed total number of crossovers Ctot = 2C(5 — 4) 
performed by the GEM algorithm of Sec. IIII Bl is not 
optimal in view of the hardness variations observed. Ad- 
ditionally, one needs to tune C for best performance. It 
turns out that, in fact, the number of crossovers can be 
determined automatically and on the run. This is done 
by comparing each pair of configurations generated by 
crossover and potential replacement of the parents: if 
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FIG. 11: (Color online) Correlation diagram of failure prob- 
abilities pf for GEM runs of length S — 64 and the total 
number of crossovers Ctot in ODGEM (overlap-driven GEM) 
runs with T> = 0.128 for 16 x 16 disorder samples. The ver- 
tical line shows the total number Ctot = 960 of crossovers for 
the GEM runs. 



they are too similar, one of them is removed from the 
population. The similarity is here again measured by the 
optimized scalar overlap q af} resulting from Eq. (flU]) , us- 
ing a cutoff g max = 1 — 2V/L 2 . For an Ising spin glass, 
T> would correspond to the number of lattice sites where 
the two configurations disagree. For the algorithm of 
Sec. IIII B[ this means that the loop over c in lines 5 and 
24 as well as the halving step in lines 25-32 are removed, 
whereas the following instructions are inserted after line 
23: 

24: if Overlap (a, (3) > 1 - 2V/L 2 then 
25: remove configuration a 

26: S <— S — 1 

27: end if 

This modified algorithm is referred to here as ODGEM 
(overlap-driven GEM). This prescription ensures max- 
imal genetic diversity at all times, and only members 
able to produce "novel" configurations in crossover are 
retained in the population. As a consequence, the to- 
tal number of crossovers is no longer fixed, but depends 
on the used disorder configuration. It appears plausible 
that hard disorder samples with many conflicting solu- 
tions close to the global minimum retain genetic diver- 
sity longer than easy samples. Figure [IT] shows a cor- 
relation plot between the failure probabilities pf in the 
original GEM and the total number of crossovers Ctot 
in the ODGEM approach (which is proportional to the 
total computational effort). The Pearson correlation co- 
efficient [Hi, 

_ ( |.Y- WH r -(>■)» 

(7 A" Oy 

is found to be p Pf ,c tot = 0.625(15), indicating clear but 
not perfect correlation. As is seen from Fig. [Til most 



FIG. 12: (Color online) Average failure probabilities for the 
GEM and ODGEM techniques on 16 x 16 samples as a func- 
tion of the total number Ctot of crossovers performed. The 
lines show fits of the functional form (|13[) to the data. 

samples with high failure rates in GEM receive an in- 
creased effort in ODGEM, and only a few cases are 
missed. The choice of the overlap cutoff V must ensure 
that when comparing configurations genetic "equality" 
is assumed only when no relative, non-trivial excitation 
exists. Values of V = 0.1 to V — 1 are found appropriate 
here. 

Since the ODGEM method detects a significant pro- 
portion of hard samples, the distribution of the total 
number C to t of crossovers acquires itself a heavy tail, re- 
flecting the hardness distribution described by Eq. (|15|) . 
This empirical distribution is found to be well modeled 
by the extreme- value shape (|T5|) with £ close to zero, i.e., 
a Gumbel form. As is evident from the data presented in 
Fig. I12[ this leads to an improved average performance of 
the ODGEM compared to the GEM technique, increas- 
ing with the total effort invested. Compared to GEM, 
ODGEM invests more effort in the hard samples and less 
in the easy ones, as appears adequate. Direct inspec- 
tion of the histogram H(pf) of failure probabilities for 
the ODGEM method reveals that, indeed, the number 
of instances with large failure probabilities are dramati- 
cally reduced as compared to the GEM data of the same 
population size presented in Fig. [5J 



2. Hardness indicators 

The possibility of deciding about sample hardness a 
priori could largely increase the efficiency and reliabil- 
ity of the GEM (and any other ground-state search) 
method. It is an open problem whether for spin-glass 
systems there are microscopic sample properties signifi- 
cantly more easily computable than the ground state it- 
self which feature a strong correlation to sample hardness 
(with respect to a given method) @,[z3,[zl- To investi- 
gate this question in the context of the GEM approach 
and the 2D bimodal XY spin glass, I analyzed corre- 
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lations between the failure probabilities pf and a large 
number of observablcs of the disorder sample, the candi- 
date ground states found, and properties of the popula- 
tion of configurations in the genetic algorithm. A num- 
ber of resulting correlation coefficients as estimated for 
16 x 16 samples from GEM and ODGEM runs with pop- 
ulation size S = 64 are collected in Tab. fl] 

Regarding easily computable properties of the disorder 
realizations at hand, one finds some moderately signifi- 
cant, positive correlation between pf in GEM runs and 
the number of frustrated plaquettes, indicating increased 
hardness for highly frustrated samples. Vice versa, an 
increasing average size of unfrustrated regions leads to 
more and more success in finding ground states. Appar- 
ently, these correlations are successfully taken into ac- 
count by the modified algorithm ODGEM, where no sig- 
nificant correlations between pf and sample properties 
are left. Also, some properties of the computed (candi- 
date) ground states correlate with success probabilities. 
In particular, a larger ground-state energy, indicating in- 
creased frustration, is accompanied by a larger number 
of failures. Also, larger than average values of the con- 
figurational r.m.s. chirality k [74j . 

K 2 L 2 =]T{ J2 signJ^x S,] z } 2 , (20) 
□« (i,3)eQ„ 

as well as the non-collinearity Q (75j . 

Q 4 i 4 /2 = E E l^iXS/, (21) 

correlate with larger failure probabilities and hence in- 
dicate harder samples for the method. Note that the 
dcfinition (|20[) is specific to the case of planar spins. For 
the Heisenberg model, the chirality is, instead, cubic in 
the spin variables [76| , Again, such correlations are not 
(highly) significant for the ODGEM technique, indicat- 
ing that the corresponding samples automatically receive 
higher computational effort there. 

The largest correlations with the failure probability are 
seen for properties of the population of spin configura- 
tions inside of the GEM or ODGEM runs. The rather 
strong inverse correlation between the average optimized 
overlap (q a/3 ) of configurations and pf in GEM runs 
shows that a homogeneous (i.e., large overlap) popula- 
tion occurs for a clear-cut, more easily accessible ground 
state. Such homogeneity also implies that larger domains 
are being identified in the BED decomposition. On the 
contrary, a large number of successful replacements of 
parents by better offspring indicates stronger competi- 
tion of candidate ground states, resulting in more fail- 
ures. These correlations related to population remain, 
although weakened, in the maximum-diversity version 
ODGEM. Consequently, devoting additional effort to dis- 
order configurations singled out by these population char- 
acteristics in ODGEM runs will additionally reduce fail- 
ure rates for hard samples, as I now discuss. 



Observable 


PGEM 


PODGEM 


No. AF bonds 


0.020(36) 


0.019(34) 


JNo. frustrated plaquettes 


0.104(33) 


—0.028(33) 


Size of unfrust. clusters 


—0.129(38) 


—0.012(37) 


Energy 


0.226(33) 


0.061(34) 


Magnetization 


0.021(34) 


0.024(36) 


Chirality 


0.137(34) 


0.001(33) 


Non-collinearity 


0.266(33) 


0.101(33) 


Plain overlap 


-0.087(31) 


0.059(32) 


Optimized overlap 


-0.414(25) 


-0.339(28) 


BED domain size 


-0.383(22) 


0.237(34) 


Parent replacements 


0.328(29) 


0.297(34) 



TABLE I: Estimated correlation coefficients between the fail- 
ure probabilities pj in GEM and ODGEM runs and various 
observables of the disorder realizations, the actual ground- 
state configurations found, and the population of configura- 
tions within the genetic algorithm. 



3. Repeated runs 

Allocating such additional effort to allegedly hard sam- 
ples typically means performing additional, statistically 
independent runs to finally pick the lowest-energy state 
found as the final result. As discussed above in Sec. lIVBl 
however, the decrease in failure probability expected 
from such a combination of several runs as described by 
Eq. (TT2]) is identical to the effect of performing a sin- 
gle computation with a larger population (and this stays 
true for the modified ODGEM technique, at least to a 
good approximation). In contrast to single, more expen- 
sive runs for all disorder realizations, however, repeated 
runs allow to treat individual realizations differently, in 
accord with the heavy-tailed distribution of sample hard- 
ness observed in Fig.H 

Even disregarding the use of the hardness indicators 
discussed in the previous Section, however, it turns out 
to be beneficial to replace runs of length T by a num- 
ber n of shorter runs of length T/n: within the range of 
validity of Eq. (fT2")) , the total probability p f not to find 
the ground state remains unchanged. For an easy sample 
with small pf, all but a small fraction of the n runs will 
end with a state of the same (namely the ground state) 
energy. For hard samples with larger pf, however, states 
with different energies will result from a sizable fraction 
of runs, even if none of them is a ground state. In other 
words, the structure of low-lying excited states typically 
results in the appearance of a variety of energies for sam- 
ples where a ground state is not found. By reacting to 
these events with additional runs for the affected samples 
until the same minimum energy has been found a certain 
number no of times, the average failure probability (pf) 
can be further decreased. 

To demonstrate the power of this extension, n — 3 runs 
of length T = S = 32 for 1 000 samples of size 16 x 16 were 
performed. For 722 samples, all three runs ended with 
the same minimum energy, which was consequently ac- 
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cepted as estimate for the ground-state energy. For the 
remaining replica, an additional run with S = 32 was 
performed, which settled 74 of the "questionable" cases. 
This scheme was iterated until for all 1 000 samples the 
state of lowest energy had been found uq = 3 times. 
The average effort for this computation corresponded to 
S w 115 (i.e., 3.6 runs of length S = 32), but the total 
number of missed ground states corresponded to that of 
runs with S « 240 (or 7.5 runs of length S = 32). This 
type of computation can, of course, be favorably com- 
bined with the correlation results of Tab. U to perform 
a certain number of additional runs for samples where 
certain observable values indicate especially large failure 
probability pf. Since the improvement effected by this 
addition depends on the structure of low-lying excited 
states of the model, it unfortunately cannot be quanti- 
fied in general. 



V. CONCLUSIONS 

I have presented a novel optimization heuristic for find- 
ing numerically exact ground states of two-dimensional 
spin-glass systems with continuous spins with high reli- 
ability. Embedding Ising spins into the continuous ro- 
tators, this exponentially hard optimization problem is 
being related to the polynomial problem of finding Ising 
ground states on planar graphs via Edmonds' algorithm 
[43j for solving minimum-weight perfect matching prob- 
lems. Due to a history dependence of effective cou- 
pling constants, this technique exhibits metastability on 
a low-energy subset of the metastable states of a single 
spin-flip zero-temperature quench. In contrast to simu- 
lated annealing and similar techniques, however, embed- 
ded matching has the crucial advantage of being strictly 
downhill in energy. To find true ground states, embed- 
ded matching is inserted as minimization procedure in 
a genetic algorithm specially tailored to the spin-glass 
ground-state problem. The essential component is here 
given by a properly chosen crossover operation exchang- 
ing automatically determined domains of rigidly locked 
spins between the parent replica, thus preserving the 
good optimization achieved at intermediate stages inside 
of domains and effectively allowing the method to di- 
rectly operate on the manifold of metastable states. 

This combination of techniques resulting in the ge- 
netic embedded matching (GEM) method outperforms 
general-purpose approaches such as simulated annealing 
by orders of magnitude: a success probability p s > 1% 
could not be achieved at all with reasonable computa- 
tional effort with simulated annealing runs for the 16 x 16 
bimodal XY samples considered for performance com- 
parison. Due to the generally strong corrections to scal- 
ing present in spin-glass systems, the extension in ac- 
cessible system sizes effected by the GEM approach over 
general-purpose techniques turns out to be crucial for the 
understanding of the as ym ptotic behavior of the spin- 
glass phase, cf. Refs. [39l.l40l|. The distribution of success 



probabilities of the GEM technique can be understood 
from the decomposition theorem (TT2"]) of failure proba- 
bilities. The thus estimated distribution over disorder 
replica of minimum required runtimes or population sizes 
is perfectly described by a Frechet distribution known 
from extreme- value theory, which is plausible given that 
sample hardness is determined by the hardest of a num- 
ber of effective two-level systems describing the energy 
landscape of a disorder realization. Due to the heavy 
tail of this distribution, the exponential divergence of 
average computational effort with system size expected 
from non-polynomial optimization problems is accom- 
panied by an increasing spread in sample hardness im- 
peding an appropriate choice of optimization parameters 
common to all disorder samples. The variant approach 
ODGEM automatically maximizing genetic diversity by 
monitoring configurational overlap, reduces the severity 
of this spread by devoting additional effort to hard sam- 
ples. Hard samples can also be detected by indicator ob- 
servables of the disorder and low-energy configurations as 
well as the population in the genetic algorithm in order 
to decrease the failure probability in these cases. The 
decomposition property (|12p finally allows to break up 
computations in smaller units which, besides allowing to 
further reduce the average effort required for a given suc- 
cess probability, makes the method ideally suitable for 
computations on parallel workstation and Beowulf clus- 
ters. 

Note that for systems with degenerate ground states, 
the GEM and ODGEM methods as non-equilibrium tech- 
niques do not yield these different states with prob- 
abilities proportional to the corresponding Boltzmann 
factors, such that in these cases an additional post- 
processing of the states found would become necessary 
[771 ] . Recent evidence suggests, however, that such de- 
generacies might be very unusual in disordered systems 
with continuous spins (33 . [loT ] . The performance analysis 
presented here focused on the bimodal XY spin glass on 
the square lattice. The method straightforwardly gen- 
eralizes to any other nearest-neighbor O(n) spin model 
(fT]) on planar graphs and to arbitrary disorder distribu- 
tions. Specifically, the case of Gaussian bond distribu- 
tion can be treated with similar efficiency. Due to the 
use of embedded matching as minimization component, 
the present from of the technique is limited to planar, 
two-dimensional lattices and the case of zero field. Other 
minimization techniques might be used in lieu of embed- 
ded matching inside of the genetic algorithm to tackle 
spin glasses in three dimensions or including magnetic 
fields. 
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